
clear;
close all;
%% select channel
P=".\data";
GAN_Path=P;
flag=1; % 1: NOF1; 2: NSC1
% channel1='Generated_data\NOF1_256';
% channel2='Generated_data\NCS1_256'; 
% filepath1 = fullfile(GAN_Path, channel1, 'mat');
% filepath2 = fullfile(GAN_Path, channel2, 'mat');
% 
% temp=1;
endloop=60;

%% 定义信道参数
%--------------------------------------------------------------------------

M_mod = 4;                         % 4QAM
k = log2(M_mod);                    % Bits/symbol
M = 512;
Mzp = 180;
B = 5000;                      % 带宽 Hz
scs = 10;                      % 子载波间隔 Hz
ofdmSym = 16;                   % ofdm信号数量 /帧 
fc = 12500;
c = 1500;                                                     % speed of light (m/s)
N = ofdmSym;

Fn = dftmtx(N);       % Generate the DFT matrix
Fn = Fn./norm(Fn);    % normalize the DFT matrix

Fm = dftmtx(M);       % Generate the DFT matrix
Fm = Fm./norm(Fm);    % normalize the DFT matrix

Wn=fwht(eye(N));      % Generate the WHT matrix
Wn=Wn./norm(Wn);      % normalize the WHT matrix

Md = M-Mzp;


%% 发送信号
data_grid = zeros(M,N);
data_grid(1:Md,:) = 1;
N_syms_perfram = sum(sum(data_grid));

delta_f = scs;
T = 1/delta_f;

             
%%                       OTFS BER Calculation

     data_na= randi([0,1],N_syms_perfram*k,1);
     data=qammod(reshape(data_na,k,Md*N), M_mod,'gray','InputType','bit');
     X = Gen_2D_data_grid(N,M,data,data_grid);

%% detect
    eng_sqrt = (M_mod==2)+(M_mod~=2)*sqrt((M_mod-1)/6*(2^2));
    SNR_dB = -6:3:15;
    SNR = 10.^(SNR_dB/10);
    noise_var_sqrt = sqrt(1./SNR);
    sigma_2 = abs(sqrt(eng_sqrt)*noise_var_sqrt).^2;

    errorRate1 = comm.ErrorRate('ResetInputPort',true);
    berOTFS1 = zeros(length(SNR),3);

    errorRate2 = comm.ErrorRate('ResetInputPort',true);
    berOFDM1 = zeros(length(SNR),3);

    errorRate3 = comm.ErrorRate('ResetInputPort',true);
    berOTSM1 = zeros(length(SNR),3);

    errorRate4 = comm.ErrorRate('ResetInputPort',true);
    berOCDM1 = zeros(length(SNR),3);

    errorRate5 = comm.ErrorRate('ResetInputPort',true);
    berAFDM1 = zeros(length(SNR),3);

    
    

    c2_AFDM = 1/(2*M); 
    for temp = 1:endloop
        disp(temp)
        channel1='Generated_data\NOF1_256';
        filepath1 = fullfile(GAN_Path, channel1, 'mat');
        filename1 = fullfile(filepath1, ['mat_' sprintf('%03d', temp) '.mat']);
        var_struct1=load(filename1);
        name_cell=fieldnames(var_struct1);
        mat1=getfield(var_struct1,char(name_cell));
        h1_real=mat1(:,:,1);
        h1_img=mat1(:,:,2);
        h1=h1_real+1i*h1_img;
        [~,G1,maxDelay,maxDoppler]=channel_acc2(h1,M,N,1);
        c1_AFDM = (2*(maxDoppler) + 1)/(2*M); % Satisfying the orthogonality condition
        % maxDoppler=20;       % Maximum normalised digital velocity
        % if (2*(maxDoppler )*(maxDelay + 1)) + (maxDelay) > M
        %     fprintf("Integer orthogonality of the AFDM is not satisfied!\n");
        % end

        for m = 1:length(SNR)
        
        
            %% modulation%%%%
        X1_tilda=X*Fn';
        s1= reshape(X1_tilda,N*M,1);

        X2_tilda=Fm'*X;
        s2= reshape(X2_tilda,N*M,1);
    
        X3_tilda=X*Wn;
        s3= reshape(X3_tilda,N*M,1);
    
        X4_tilda=OCDMmod(X,M);
        s4= reshape(X4_tilda,N*M,1);
    
        X5_tilda=AFDMmod(X,c1_AFDM,c2_AFDM);
        s5= reshape(X5_tilda,N*M,1);

        noise= sqrt(sigma_2(m)/2)*(randn(size(s1)) + 1i*randn(size(s1)));

        r1=G1*s1;
        r1_n=r1+noise;

        r2=G1*s2;
        r2_n=r2+noise;

        r3=G1*s3;
        r3_n=r3+noise;

        r4=G1*s4;
        r4_n=r4+noise;

        r5=G1*s5;
        r5_n=r5+noise;


        [est_LMMSE_r1_n,~] = LMMSE_detector_OTFS(N,M,M_mod,sigma_2(m),data_grid,r1_n,G1); 
        [est_LMMSE_r2_n,~] = LMMSE_detector_OFDM(N,M,M_mod,sigma_2(m),data_grid,r2_n,G1); 
        [est_LMMSE_r3_n,~] = LMMSE_detector_OTSM(N,M,M_mod,sigma_2(m),data_grid,r3_n,G1); 
        [est_LMMSE_r4_n,~] = LMMSE_detector_OCDM(N,M,M_mod,sigma_2(m),data_grid,r4_n,G1); 
        [est_LMMSE_r5_n,~] = LMMSE_detector_AFDM(N,M,M_mod,sigma_2(m),data_grid,r5_n,G1,c1_AFDM,c2_AFDM);
        
        berOTFS1(m,:)=berOTFS1(m,:)+errorRate1(data_na,est_LMMSE_r1_n,1).';
        berOFDM1(m,:)=berOFDM1(m,:)+errorRate2(data_na,est_LMMSE_r2_n,1).';
        berOTSM1(m,:)=berOTSM1(m,:)+errorRate3(data_na,est_LMMSE_r3_n,1).';
        berOCDM1(m,:)=berOCDM1(m,:)+errorRate4(data_na,est_LMMSE_r4_n,1).';
        berAFDM1(m,:)=berAFDM1(m,:)+errorRate5(data_na,est_LMMSE_r5_n,1).';
            
        end
      
    end

%% figure
% berOTFS1(~berOTFS1)=1e-12;
% berOTFS2(~berOTFS2)=1e-12;


figure
berOFDM1(:,1)=berOFDM1(:,2)./berOFDM1(:,3);
berOCDM1(:,1)=berOCDM1(:,2)./berOCDM1(:,3);
berAFDM1(:,1)=berAFDM1(:,2)./berAFDM1(:,3);
berOTSM1(:,1)=berOTSM1(:,2)./berOTSM1(:,3);
berOTFS1(:,1)=berOTFS1(:,2)./berOTFS1(:,3);

save('./ber_data/main_channel_avg_NOF1_gen_ber.mat','berOFDM1','berOCDM1','berAFDM1','berOTSM1','berOTFS1');

% berOTFS1(~berOTFS1)=1e-6;
% berOFDM1(~berOFDM1)=1e-6;
% berOTSM1(~berOTSM1)=1e-6;
% berOCDM1(~berOCDM1)=1e-6;
% berAFDM1(~berAFDM1)=1e-6;


semilogy(SNR_dB,berOFDM1(:,1),'Marker','o','LineWidth',2,'Color',[0 0 1]);
hold on;
semilogy(SNR_dB,berOCDM1(:,1),'Marker','+','LineWidth',2,'Color',[0,0,0]);
semilogy(SNR_dB,berAFDM1(:,1),'Marker','*','LineWidth',2,'Color',[0.8500 0.3250 0.0980]);
semilogy(SNR_dB,berOTSM1(:,1),'Marker','x','LineWidth',2,'Color',[0 1 0]);
semilogy(SNR_dB,berOTFS1(:,1),'Marker','square','LineWidth',2,'Color',[1 0 0]);   

xlim([-6,15]);

ylabel('BER');
xlabel('SNR/dB');
legend( 'OFDM',"OCDM","AFDM","OTSM","OTFS");
title("NOF1生成信道BER compare")
grid on;
hold off;
